1. laboratorijska vježba

Multivarijatna analiza podataka

Verzija: 1.0

1. Uvod i upute za predaju

Cilj ove laboratorijske vježbe je primijeniti osnovne koncepte multivarijatne analize podataka, istražiti podatke te ispitati hipoteze. Preduvjet za rješavanje vježbe je osnovno znanje programskog jezika R i rad s R Markdown dokumentima. Sama vježba je koncipirana kao projekt u kojem istražujete i eksperimentirate koristeći dane podatke - ne postoji nužno samo jedan točan način rješavanja svakog podzadatka.

Rješavanje vježbe svodi se na čitanje uputa u tekstu ovog dokumenta, nadopunjavanje blokova kôda (možete dodavati i dodatne blokove kôda ukoliko je potrebno) i ispisivanje rezultata (u vidu ispisa iz funkcija, tablica i grafova). Vježbu radite samostalno, a svoje rješenje branite na terminima koji su vam dodijeljeni u kalendaru. Pritom morate razumjeti teorijske osnove u okviru onoga što je obrađeno na predavanjima i morate pokazati da razumijete sav kôd koji ste napisali.

Vaše rješenje potrebno je predati u sustav Moodle u obliku dvije datoteke:

1. Ovaj .Rmd dokument s Vašim rješenjem (naziva IME_PREZIME_JMBAG.rmd),
   
2. PDF ili HTML dokument kao izvještaj generiran iz vašeg .Rmd rješenja (također naziva IME_PREZIME_JMBAG).

Rok za predaju je 3. travnja 2022. u 23:59h. Podsjećamo da bodovi iz laboratorijskih vježbi ulaze i u bodove na ispitnom roku, te da je za polaganje predmeta potrebno imati barem 50% ukupnih bodova iz laboratorijskih vježbi. Nadoknade laboratorijskih vježbi neće biti organizirane. Za sva dodatna pitanja svakako se javite na email adresu predmeta: map@fer.hr.

2. Podatkovni skup

Podatkovni skup koji će biti razmatran u vježbi sadrži bodove studenata na jednom fakultetskom kolegiju. Svakom studentu upisani su bodovi iz dviju laboratorijskih vježbi (LAB), pet zadataka međuispita (MI), pet zadataka završnog ispita (ZI), pet zadataka ispitnog roka (IR) i kojoj grupi predavanja pripadaju (Grupa).

Studenti mogu položiti kolegij kontinuiranim putem ili na ispitnom roku. Kontinuirani put sastoji se od bodova s laboratorijskih vježbi, međuispita i završnog ispita. Kronološki, 1. laboratorijska vježba održana je prije međuispita, dok je 2. laboratorijska vježba održana između međuispita i završnog ispita. Ispitni rok održan je nakon završnog ispita. Ako student polaže predmet na ispitnom roku, gledaju se samo bodovi s ispitnog roka. Ukupan broj bodova je 100, a bodovi su raspodijeljeni na sljedeći način:

• Kontinuirana nastava:
  • LAB: 20 bodova (0-10 svaka vježba)
  • MI : 40 bodova (0-8 svaki zadatak)
  • ZI : 40 bodova (0-8 svaki zadatak)
• Ispitni rok:
  • IR : 100 bodova (0-20 svaki zadatak)

Za prolazak kolegija potrebno je skupiti više od 50 bodova i izaći na obje laboratorijske vježbe (izlazak na vježbe nužan je uvjet i za polaganje ispitnog roka, iako se bodovi ne prenose). Ako student nije pristupio pripadajućem ispitu/laboratorijskoj vježbi, nije upisan podatak (što nije isto kao i 0 bodova).

3. Priprema podataka i eksploratorna analiza

U ovom dijelu vježbe potrebno je učitati podatke i napraviti osnovnu eksploratornu analizu podataka.

3.1 Učitavanje podataka

Učitajte podatkovni skup iz datoteke studenti.csv i pripremite podatke za analizu. Pritom obratite pozornost na sljedeće:

• Provjerite jesu li sve varijable očekivanog tipa,
• Provjerite jesu li vrijednosti unutar zadanog raspona (s obzirom na gore opisano bodovanje),
• Provjerite zadovoljavaju li bodovi gore opisane uvjete predmeta,
• Za nedostajuće podatke ispitajte jesu li opravdani te odaberite i primijenite tehniku upravljanja nedostajućim podatcima.

Nakon što su podatci pripremljeni, analizirajte i ispišite deksriptivne statistike varijabli.

grades <- read.csv('studenti.csv')
summary(grades)

##       MI_1            MI_2             MI_3            MI_4      
##  Min.   :4.000   Min.   : 0.000   Min.   :0.000   Min.   :0.500  
##  1st Qu.:6.500   1st Qu.: 4.500   1st Qu.:3.500   1st Qu.:3.500  
##  Median :7.000   Median : 6.000   Median :5.000   Median :4.000  
##  Mean   :6.912   Mean   : 5.847   Mean   :4.926   Mean   :4.007  
##  3rd Qu.:7.500   3rd Qu.: 7.500   3rd Qu.:6.500   3rd Qu.:4.500  
##  Max.   :8.000   Max.   :18.000   Max.   :8.000   Max.   :7.000  
##                                                                  
##      MI_5              LAB_1                ZI_1             ZI_2      
##  Length:500         Length:500         Min.   :-3.000   Min.   :3.000  
##  Class :character   Class :character   1st Qu.: 4.500   1st Qu.:5.500  
##  Mode  :character   Mode  :character   Median : 6.000   Median :6.000  
##                                        Mean   : 5.811   Mean   :5.985  
##                                        3rd Qu.: 7.500   3rd Qu.:6.500  
##                                        Max.   : 8.000   Max.   :8.000  
##                                                                        
##       ZI_3            ZI_4           ZI_5              LAB_2          
##  Min.   :0.000   Min.   :0.000   Length:500         Length:500        
##  1st Qu.:2.500   1st Qu.:2.500   Class :character   Class :character  
##  Median :4.000   Median :3.000   Mode  :character   Mode  :character  
##  Mean   :4.009   Mean   :2.997                                        
##  3rd Qu.:5.500   3rd Qu.:3.500                                        
##  Max.   :8.000   Max.   :5.500                                        
##                                                                       
##       IR_1            IR_2            IR_3            IR_4      
##  Min.   : 0.00   Min.   : 0.00   Min.   : 0.00   Min.   : 0.00  
##  1st Qu.:13.50   1st Qu.:12.50   1st Qu.:13.50   1st Qu.: 8.00  
##  Median :15.00   Median :14.50   Median :14.50   Median :11.25  
##  Mean   :15.25   Mean   :14.07   Mean   :14.34   Mean   :11.14  
##  3rd Qu.:17.62   3rd Qu.:16.00   3rd Qu.:15.50   3rd Qu.:14.12  
##  Max.   :20.00   Max.   :20.00   Max.   :18.50   Max.   :20.00  
##  NA's   :400     NA's   :400     NA's   :400     NA's   :400    
##       IR_5            Grupa      
##  Min.   : 0.000   Min.   :1.000  
##  1st Qu.: 5.000   1st Qu.:1.000  
##  Median : 6.500   Median :2.000  
##  Mean   : 6.305   Mean   :2.006  
##  3rd Qu.: 7.500   3rd Qu.:3.000  
##  Max.   :11.500   Max.   :3.000  
##  NA's   :400

Vidimo da MI_5, LAB_1, ZI_5, i LAB_2 trebamo prebaciti u numerički tip podatka.

grades[ , colnames(grades)] <- apply(grades[ , colnames(grades)], 2, # Specify own function within apply
                                     function(x) as.numeric(as.character(x)))

## Warning in FUN(newX[, i], ...): NAs introduced by coercion

## Warning in FUN(newX[, i], ...): NAs introduced by coercion

## Warning in FUN(newX[, i], ...): NAs introduced by coercion

## Warning in FUN(newX[, i], ...): NAs introduced by coercion

head(grades)

sapply(grades, class)

##      MI_1      MI_2      MI_3      MI_4      MI_5     LAB_1      ZI_1      ZI_2 
## "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" 
##      ZI_3      ZI_4      ZI_5     LAB_2      IR_1      IR_2      IR_3      IR_4 
## "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" 
##      IR_5     Grupa 
## "numeric" "numeric"

summary(grades)

##       MI_1            MI_2             MI_3            MI_4      
##  Min.   :4.000   Min.   : 0.000   Min.   :0.000   Min.   :0.500  
##  1st Qu.:6.500   1st Qu.: 4.500   1st Qu.:3.500   1st Qu.:3.500  
##  Median :7.000   Median : 6.000   Median :5.000   Median :4.000  
##  Mean   :6.912   Mean   : 5.847   Mean   :4.926   Mean   :4.007  
##  3rd Qu.:7.500   3rd Qu.: 7.500   3rd Qu.:6.500   3rd Qu.:4.500  
##  Max.   :8.000   Max.   :18.000   Max.   :8.000   Max.   :7.000  
##                                                                  
##       MI_5          LAB_1            ZI_1             ZI_2      
##  Min.   :0.00   Min.   :4.000   Min.   :-3.000   Min.   :3.000  
##  1st Qu.:1.50   1st Qu.:6.500   1st Qu.: 4.500   1st Qu.:5.500  
##  Median :3.00   Median :7.000   Median : 6.000   Median :6.000  
##  Mean   :3.05   Mean   :6.995   Mean   : 5.811   Mean   :5.985  
##  3rd Qu.:4.50   3rd Qu.:7.500   3rd Qu.: 7.500   3rd Qu.:6.500  
##  Max.   :8.00   Max.   :9.500   Max.   : 8.000   Max.   :8.000  
##  NA's   :1      NA's   :2                                       
##       ZI_3            ZI_4            ZI_5           LAB_2      
##  Min.   :0.000   Min.   :0.000   Min.   :0.000   Min.   :0.500  
##  1st Qu.:2.500   1st Qu.:2.500   1st Qu.:1.500   1st Qu.:2.500  
##  Median :4.000   Median :3.000   Median :2.000   Median :3.000  
##  Mean   :4.009   Mean   :2.997   Mean   :2.019   Mean   :3.003  
##  3rd Qu.:5.500   3rd Qu.:3.500   3rd Qu.:2.500   3rd Qu.:3.500  
##  Max.   :8.000   Max.   :5.500   Max.   :5.500   Max.   :6.000  
##                                  NA's   :1       NA's   :2      
##       IR_1            IR_2            IR_3            IR_4      
##  Min.   : 0.00   Min.   : 0.00   Min.   : 0.00   Min.   : 0.00  
##  1st Qu.:13.50   1st Qu.:12.50   1st Qu.:13.50   1st Qu.: 8.00  
##  Median :15.00   Median :14.50   Median :14.50   Median :11.25  
##  Mean   :15.25   Mean   :14.07   Mean   :14.34   Mean   :11.14  
##  3rd Qu.:17.62   3rd Qu.:16.00   3rd Qu.:15.50   3rd Qu.:14.12  
##  Max.   :20.00   Max.   :20.00   Max.   :18.50   Max.   :20.00  
##  NA's   :400     NA's   :400     NA's   :400     NA's   :400    
##       IR_5            Grupa      
##  Min.   : 0.000   Min.   :1.000  
##  1st Qu.: 5.000   1st Qu.:1.000  
##  Median : 6.500   Median :2.000  
##  Mean   : 6.305   Mean   :2.006  
##  3rd Qu.: 7.500   3rd Qu.:3.000  
##  Max.   :11.500   Max.   :3.000  
##  NA's   :400

dim(grades)

## [1] 500  18

Rasponi vrijednosti MI_2, ZI_1 nisu u skladu s opisom domene.

Uklanjanje redaka s vrijednostima izvan dopuštenog raspona

grades[grades$MI_2 > 8 | grades$ZI_1 < 0,]

grades <- grades[!(grades$MI_2 > 8 | grades$ZI_1 < 0),]

Vidimo da su samo dva reda problematična pa smo ta dva reda izbacili iz dataseta.

Izbacivanje redaka s problematičnim nedostajućim vrijednostima

grades[is.na(grades$MI_5) | is.na(grades$ZI_5), ]

mi_or_zi <- colnames(grades)[colnames(grades) %like% '[M|Z]I']
grades <- grades[!rowSums(is.na(grades[, mi_or_zi])) > 0,]
dim(grades)

## [1] 496  18

Također vidimo da su svi studenti pristupili MI i ZI, no u dva retka za neki od zadataka nedostaje vrijednost. Te retke ćemo izbaciti jer je za njih broj bodova izgubljen.

Ispitivanje nedostajućih vrijednosti IR

missing_ir <- grades[rowSums(is.na(grades[,colnames(grades)[colnames(grades) %like% 'IR']])) > 0, ]
dim(missing_ir)

## [1] 399  18

Iz ove provjere vidimo da svaki student koji ima ‘NA’ iz nekog zadatka na IR ima ‘NA’ na svim zadatcima s IR, to jest da nije prisustvovao IR. Dakle, nedostajuće vrijednosti laboratorijskih vježbi i Ispitnih rokova možemo ostaviti nepromijenjene.

3.2 Korelacijska analiza

Razmotrimo studente koji su predmet položili kontinuirano. Izračunajte i vizualizirajte matricu korelacije za njihove bodove na nastavnim aktivnostima. Ponovite isto za studente koji su izašli na ispitni rok. Razmislite o zavisnosti različitih nastavnih aktivnosti koje vidite iz ovih korelacijskih matrica.

has_both_labs <- grades[!rowSums(is.na(grades[, c('LAB_1', 'LAB_2')])) > 0, ] # Izbacimo sve studente koji nisu bili na svim labosa
dim(has_both_labs) # Izbacili smo 3 studenta

## [1] 493  18

Mi_or_Zi_or_Lab <- colnames(grades)[colnames(grades) %like% '[M|Z]I|LAB']
kont_pass <- has_both_labs[rowSums(has_both_labs[, Mi_or_Zi_or_Lab]) > 50, ]
dim(kont_pass)

## [1] 399  18

cols_except_IR <- colnames(grades)[!colnames(grades) %like% 'IR']
kont_pass <- kont_pass[, cols_except_IR]
corr_matrix <- cor(kont_pass)
corr_matrix

##              MI_1        MI_2         MI_3         MI_4        MI_5       LAB_1
## MI_1   1.00000000  0.36540468  0.047202152  0.209778250 -0.03042302  0.37712960
## MI_2   0.36540468  1.00000000  0.106896046  0.012687469 -0.02014711  0.04744858
## MI_3   0.04720215  0.10689605  1.000000000  0.005550264  0.03583983  0.09525004
## MI_4   0.20977825  0.01268747  0.005550264  1.000000000  0.10812560  0.05335674
## MI_5  -0.03042302 -0.02014711  0.035839829  0.108125603  1.00000000 -0.06807928
## LAB_1  0.37712960  0.04744858  0.095250035  0.053356738 -0.06807928  1.00000000
## ZI_1  -0.08943602 -0.11397463  0.132026817 -0.070774859 -0.11250821 -0.06606030
## ZI_2  -0.07060226 -0.02383213 -0.110486873 -0.056920820  0.02684550 -0.08190418
## ZI_3  -0.11441076 -0.15955231 -0.100154730  0.018164160 -0.14861858 -0.11737454
## ZI_4  -0.11987925 -0.10158006 -0.048150360 -0.061583381 -0.03816174 -0.07561181
## ZI_5  -0.04986475 -0.02530377 -0.054419970 -0.011229705 -0.05954402 -0.06344536
## LAB_2  0.15087328 -0.01180414 -0.141086262 -0.042151271 -0.11032307  0.28457155
## Grupa -0.12836154 -0.10304016 -0.136686153 -0.031954732 -0.09749442 -0.15466617
##              ZI_1         ZI_2        ZI_3        ZI_4         ZI_5       LAB_2
## MI_1  -0.08943602 -0.070602256 -0.11441076 -0.11987925 -0.049864750  0.15087328
## MI_2  -0.11397463 -0.023832133 -0.15955231 -0.10158006 -0.025303773 -0.01180414
## MI_3   0.13202682 -0.110486873 -0.10015473 -0.04815036 -0.054419970 -0.14108626
## MI_4  -0.07077486 -0.056920820  0.01816416 -0.06158338 -0.011229705 -0.04215127
## MI_5  -0.11250821  0.026845502 -0.14861858 -0.03816174 -0.059544020 -0.11032307
## LAB_1 -0.06606030 -0.081904177 -0.11737454 -0.07561181 -0.063445358  0.28457155
## ZI_1   1.00000000  0.319867127  0.07796473  0.13516805 -0.027332048  0.41264570
## ZI_2   0.31986713  1.000000000  0.21866566  0.16015212  0.007794613  0.08676891
## ZI_3   0.07796473  0.218665662  1.00000000  0.11831387  0.031992222 -0.07700355
## ZI_4   0.13516805  0.160152121  0.11831387  1.00000000  0.057387224 -0.04719475
## ZI_5  -0.02733205  0.007794613  0.03199222  0.05738722  1.000000000 -0.07683370
## LAB_2  0.41264570  0.086768908 -0.07700355 -0.04719475 -0.076833700  1.00000000
## Grupa -0.08196034 -0.041414881 -0.08672903 -0.11554507 -0.044847212 -0.09807580
##             Grupa
## MI_1  -0.12836154
## MI_2  -0.10304016
## MI_3  -0.13668615
## MI_4  -0.03195473
## MI_5  -0.09749442
## LAB_1 -0.15466617
## ZI_1  -0.08196034
## ZI_2  -0.04141488
## ZI_3  -0.08672903
## ZI_4  -0.11554507
## ZI_5  -0.04484721
## LAB_2 -0.09807580
## Grupa  1.00000000

which(corr_matrix > 0.3 & corr_matrix != 1, arr.ind = TRUE)

##       row col
## MI_2    2   1
## LAB_1   6   1
## MI_1    1   2
## MI_1    1   6
## ZI_2    8   7
## LAB_2  12   7
## ZI_1    7   8
## ZI_1    7  12

corr_matrix[which(corr_matrix > 0.3 & corr_matrix != 1, arr.ind = TRUE)]

## [1] 0.3654047 0.3771296 0.3654047 0.3771296 0.3198671 0.4126457 0.3198671
## [8] 0.4126457

which(corr_matrix < -0.15, arr.ind = TRUE)

##       row col
## ZI_3    9   2
## Grupa  13   6
## MI_2    2   9
## LAB_1   6  13

corr_matrix[which(corr_matrix < -0.15, arr.ind = TRUE)]

## [1] -0.1595523 -0.1546662 -0.1595523 -0.1546662

Ako pogledamo korelacije između svakih pojedinih zadataka i laboratorijskih vježbi nećemo pronaći jake zavisnosti između nastavnih aktivnosti. Najjača korelacija jest 0.413 između varijabli ZI_1 te LAB_2. Moguće je da je prvi zadatak završnog ispita provjeravao znanje iz druge laboratorijske vježbe pa su studenti koji su bolje riješili taj labos također postigli bolji rezultat na tom zadatku. Pogledajmo još i korelacije ukupnog rezultata na MI i ZI te vježbi.

just_MI <- colnames(grades)[colnames(grades) %like% 'MI']
just_ZI <- colnames(grades)[colnames(grades) %like% 'ZI']
kont_pass$MI_TOTAL <- rowSums(kont_pass[, just_MI])
kont_pass$ZI_TOTAL <- rowSums(kont_pass[, just_ZI])
kont_pass_reduced <- kont_pass[, c('MI_TOTAL', 'ZI_TOTAL', 'LAB_1', 'LAB_2')]
cor(kont_pass_reduced)

##            MI_TOTAL   ZI_TOTAL      LAB_1      LAB_2
## MI_TOTAL  1.0000000 -0.2338468  0.1282353 -0.1081578
## ZI_TOTAL -0.2338468  1.0000000 -0.1534233  0.1310967
## LAB_1     0.1282353 -0.1534233  1.0000000  0.2845716
## LAB_2    -0.1081578  0.1310967  0.2845716  1.0000000

Kada gledamo ukupne rezultate korelacije su još slabije.

Provedimo isti postupak nad studentima s studentskog roka:

has_both_labs <- grades[!rowSums(is.na(grades[, c('LAB_1', 'LAB_2')])) > 0, ] # Izbacimo sve studente koji nisu bili na svim labosa
dim(has_both_labs) # Izbacili smo 3 studenta

## [1] 493  18

just_IR <- colnames(grades)[colnames(grades) %like% 'IR']
both_labs_and_IR <- has_both_labs[!(rowSums(is.na(has_both_labs[, just_IR])) > 0), ]
IR_pass <- both_labs_and_IR[rowSums(both_labs_and_IR[, just_IR]) > 50, ]
dim(IR_pass)

## [1] 88 18

corr_matrix_IR <- cor(IR_pass)
corr_matrix_IR

##                MI_1         MI_2         MI_3        MI_4        MI_5
## MI_1   1.0000000000  0.391094237  0.006099483  0.04448703 -0.02496981
## MI_2   0.3910942367  1.000000000 -0.162554392  0.07694659 -0.11758017
## MI_3   0.0060994827 -0.162554392  1.000000000 -0.03272093 -0.22641599
## MI_4   0.0444870296  0.076946592 -0.032720935  1.00000000  0.10348519
## MI_5  -0.0249698133 -0.117580174 -0.226415992  0.10348519  1.00000000
## LAB_1  0.4089645605  0.040469642 -0.109397941 -0.09140253  0.06036877
## ZI_1  -0.3489449124 -0.376308137  0.137416218 -0.27150881 -0.15030582
## ZI_2  -0.1019557677 -0.250368012  0.018457995 -0.06001195 -0.03400316
## ZI_3  -0.3445820261 -0.033454578 -0.301652033  0.04323348  0.07872828
## ZI_4   0.0274361071 -0.039345514 -0.061403955 -0.04493153 -0.10015789
## ZI_5  -0.1187774456 -0.152947637 -0.033933489 -0.10277014  0.05718938
## LAB_2 -0.0333682329 -0.242664054 -0.039736089 -0.12744759  0.06199763
## IR_1   0.4688849465 -0.039646027  0.279877733  0.08143881 -0.18866381
## IR_2  -0.1494015553  0.009758445 -0.100517702  0.36182358  0.21703880
## IR_3   0.1125968460  0.049840100 -0.256493271  0.01327931  0.03171626
## IR_4  -0.0177226975 -0.194728637  0.118016315  0.16364502  0.02119351
## IR_5  -0.0238973752  0.069431584 -0.053732918  0.02560355  0.05627318
## Grupa -0.0004252649 -0.166246845  0.012723820 -0.08555069 -0.22402451
##              LAB_1        ZI_1        ZI_2         ZI_3        ZI_4
## MI_1   0.408964561 -0.34894491 -0.10195577 -0.344582026  0.02743611
## MI_2   0.040469642 -0.37630814 -0.25036801 -0.033454578 -0.03934551
## MI_3  -0.109397941  0.13741622  0.01845799 -0.301652033 -0.06140395
## MI_4  -0.091402527 -0.27150881 -0.06001195  0.043233477 -0.04493153
## MI_5   0.060368770 -0.15030582 -0.03400316  0.078728277 -0.10015789
## LAB_1  1.000000000 -0.40368156 -0.21418178 -0.193899208  0.07113723
## ZI_1  -0.403681556  1.00000000  0.27290163  0.264866950 -0.01478738
## ZI_2  -0.214181782  0.27290163  1.00000000  0.237504099  0.30463208
## ZI_3  -0.193899208  0.26486695  0.23750410  1.000000000 -0.02108238
## ZI_4   0.071137229 -0.01478738  0.30463208 -0.021082384  1.00000000
## ZI_5   0.001096857  0.11774203 -0.02790677 -0.080100662 -0.10213531
## LAB_2  0.110189996  0.31834892 -0.01160661  0.033205240 -0.04962401
## IR_1  -0.085626143 -0.10846904 -0.04079295 -0.215695732  0.02873990
## IR_2   0.037006594 -0.09470889 -0.07888448  0.047727895 -0.06418097
## IR_3   0.219190706  0.18801488 -0.30103685  0.076291180 -0.16597663
## IR_4  -0.042083578  0.02046877  0.34653553  0.117263457 -0.04289348
## IR_5  -0.047189647 -0.06912521  0.05867845  0.008698339  0.16421871
## Grupa -0.070628071  0.05844784  0.03941188 -0.095528971 -0.01959692
##               ZI_5       LAB_2         IR_1         IR_2         IR_3
## MI_1  -0.118777446 -0.03336823  0.468884947 -0.149401555  0.112596846
## MI_2  -0.152947637 -0.24266405 -0.039646027  0.009758445  0.049840100
## MI_3  -0.033933489 -0.03973609  0.279877733 -0.100517702 -0.256493271
## MI_4  -0.102770143 -0.12744759  0.081438810  0.361823580  0.013279311
## MI_5   0.057189378  0.06199763 -0.188663807  0.217038802  0.031716262
## LAB_1  0.001096857  0.11019000 -0.085626143  0.037006594  0.219190706
## ZI_1   0.117742030  0.31834892 -0.108469044 -0.094708895  0.188014878
## ZI_2  -0.027906774 -0.01160661 -0.040792948 -0.078884481 -0.301036847
## ZI_3  -0.080100662  0.03320524 -0.215695732  0.047727895  0.076291180
## ZI_4  -0.102135314 -0.04962401  0.028739898 -0.064180971 -0.165976633
## ZI_5   1.000000000  0.02209739 -0.033281260  0.111776561  0.060243188
## LAB_2  0.022097389  1.00000000 -0.042991909  0.035849793 -0.179684022
## IR_1  -0.033281260 -0.04299191  1.000000000  0.098236973  0.050418494
## IR_2   0.111776561  0.03584979  0.098236973  1.000000000  0.005378359
## IR_3   0.060243188 -0.17968402  0.050418494  0.005378359  1.000000000
## IR_4  -0.077961398  0.08570545  0.051583465  0.106099044  0.087874626
## IR_5   0.508570489 -0.14833718 -0.008335806  0.130938454  0.027192675
## Grupa  0.091305084  0.01432996  0.233686386  0.022368290 -0.003770292
##              IR_4         IR_5         Grupa
## MI_1  -0.01772270 -0.023897375 -0.0004252649
## MI_2  -0.19472864  0.069431584 -0.1662468450
## MI_3   0.11801631 -0.053732918  0.0127238203
## MI_4   0.16364502  0.025603551 -0.0855506886
## MI_5   0.02119351  0.056273176 -0.2240245126
## LAB_1 -0.04208358 -0.047189647 -0.0706280711
## ZI_1   0.02046877 -0.069125206  0.0584478355
## ZI_2   0.34653553  0.058678445  0.0394118803
## ZI_3   0.11726346  0.008698339 -0.0955289713
## ZI_4  -0.04289348  0.164218707 -0.0195969236
## ZI_5  -0.07796140  0.508570489  0.0913050842
## LAB_2  0.08570545 -0.148337183  0.0143299582
## IR_1   0.05158346 -0.008335806  0.2336863863
## IR_2   0.10609904  0.130938454  0.0223682904
## IR_3   0.08787463  0.027192675 -0.0037702916
## IR_4   1.00000000 -0.141255989  0.0455580520
## IR_5  -0.14125599  1.000000000  0.0714675992
## Grupa  0.04555805  0.071467599  1.0000000000

which(corr_matrix_IR > 0.4 & corr_matrix_IR != 1, arr.ind = TRUE)

##       row col
## LAB_1   6   1
## IR_1   13   1
## MI_1    1   6
## IR_5   17  11
## MI_1    1  13
## ZI_5   11  17

corr_matrix_IR[which(corr_matrix_IR > 0.4 & corr_matrix_IR != 1, arr.ind = TRUE)]

## [1] 0.4089646 0.4688849 0.4089646 0.5085705 0.4688849 0.5085705

which(corr_matrix_IR < -0.3, arr.ind = TRUE)

##       row col
## ZI_1    7   1
## ZI_3    9   1
## ZI_1    7   2
## ZI_3    9   3
## ZI_1    7   6
## MI_1    1   7
## MI_2    2   7
## LAB_1   6   7
## IR_3   15   8
## MI_1    1   9
## MI_3    3   9
## ZI_2    8  15

corr_matrix_IR[which(corr_matrix_IR < -0.3, arr.ind = TRUE)]

##  [1] -0.3489449 -0.3445820 -0.3763081 -0.3016520 -0.4036816 -0.3489449
##  [7] -0.3763081 -0.4036816 -0.3010368 -0.3445820 -0.3016520 -0.3010368

Najjača korelacija koju zamjećujemo jest 0.501 između IR_5 te ZI_5, to možda ukazuje da ti zadatci ispituju istu kompetenciju.

just_MI <- colnames(grades)[colnames(grades) %like% 'MI']
just_ZI <- colnames(grades)[colnames(grades) %like% 'ZI']
IR_pass$MI_TOTAL <- rowSums(IR_pass[, just_MI])
IR_pass$ZI_TOTAL <- rowSums(IR_pass[, just_ZI])
IR_pass$IR_TOTAL <- rowSums(IR_pass[, just_IR])
kont_pass_reduced <- IR_pass[, c('MI_TOTAL', 'ZI_TOTAL', 'IR_TOTAL', 'LAB_1', 'LAB_2')]
cor(kont_pass_reduced)

##             MI_TOTAL   ZI_TOTAL     IR_TOTAL        LAB_1       LAB_2
## MI_TOTAL  1.00000000 -0.4471485  0.131488736  0.070370288 -0.19070591
## ZI_TOTAL -0.44714846  1.0000000  0.047675800 -0.333717672  0.16614769
## IR_TOTAL  0.13148874  0.0476758  1.000000000 -0.004886852 -0.04615177
## LAB_1     0.07037029 -0.3337177 -0.004886852  1.000000000  0.11019000
## LAB_2    -0.19070591  0.1661477 -0.046151774  0.110189996  1.00000000

Kada gledamo ukupne brojeve bodova ne primjećujemo jake korelacije. Ipak, korelacija MI_TOTAL i ZI_TOTAL iznosi -0.447 što može biti uzrokovano time da studenti koji polože predmet na roku “pobace” na jednom od kolokvija, (da su ostvarili dobar rezultat na oba kolokvija ne bi izašli na rok).

Prikažite upareni graf za zadatke s ispitnog roka. Na dijagonalama prikažite empirijsku distribuciju podataka, a na elementima izvan dijagonala prikažite grafove raspršenja za parove varijabli. Razmislite o karakteristikama grafova i razmislite postoje li primjeri koji odskaču od ostalih.

IR_data <- grades[!(rowSums(is.na(grades[, just_IR])) > 0), ]
IR_data <- IR_data[, just_IR]
ggpairs(IR_data, diag = list(continuous = "barDiag"))+theme_bw()+ geom_histogram(aes(y = ..density..))

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Vidimo da prva tri zadatka s roka imaju puno uže distribucije od preostalih zadataka IR_4 i IR_5 koji poprimaju dosta šire raspone vrijednosti. Također uočavamo da neke distribucije poprimaju ne zvonolike oblike: IR_1, te bimodalna distribucija IR_4. Vrlo je vjerojatno da točka u donjem lijevom kutu svakog grafa pripada istom retku, to ćemo provjeriti kasnije.

3.3 Statistička udaljenost

Izračunajte procjene vektora očekivanja i matrice kovarijance za zadatke s ispitnog roka, kao i statističke udaljenosti svih primjera u odnosu na procijenjeno očekivanje i kovarijancu. Ispitajte postoje li stršeće vrijednosti koje su statistički značajne.

colMeans(IR_data)

##      IR_1      IR_2      IR_3      IR_4      IR_5 
## 15.335052 14.128866 14.376289 11.072165  6.304124

cov(IR_data)

##           IR_1     IR_2      IR_3       IR_4       IR_5
## IR_1 11.285009 3.649082 3.0601106  2.9026525  0.7642290
## IR_2  3.649082 8.928533 2.3338166  3.9515410  1.2859214
## IR_3  3.060111 2.333817 5.5314111  2.3319373  0.9807238
## IR_4  2.902652 3.951541 2.3319373 21.8332796 -0.4596757
## IR_5  0.764229 1.285921 0.9807238 -0.4596757  4.2268578

stat_dist <- mahalanobis(IR_data, colMeans(IR_data), cov(IR_data))
boxplot(stat_dist, main="Mahalanobis distances of IR_data", ylab="Mahalanobis")

hist(stat_dist)

IR_data[stat_dist > 20,]

stat_dist[stat_dist > 20]

##       35 
## 50.77283

Vidimo da postoji samo jedan outlier po statističkoj udaljenosti te je to upravo redak kojeg smo detektirali i u prošlom zadatku. Ipak, ovaj redak vjerojatno nije kriv, već samo predstavlja studenta koji je ostavio sve zadatke na roku praznima. Ovakav događaj je sasvim moguć pa taj redak nećemo izbaciti iz analize, ali ćemo ga imati na umu pri provedbi daljnjih istraživanja.

4. Analiza podataka

4.1 Vizualizacija i deskriptivna statistika

Analizirajte u podatcima sljedeća istraživačka pitanja, koristeći odgovarajuće vizualizacije i deskriptivne statistike ili druge tehnike (dodatno možete provesti i statistički test - nije obavezno).

• Imaju li grupe utjecaj na ukupne bodove iz kontinuirane nastave (postoje li grupe koje su uspješnije od ostalih)? Vrijedi li isto za bodove na roku?

grades_no_na <- grades[!(rowSums(is.na(grades[ ,Mi_or_Zi_or_Lab])) > 0), ]
grades_no_na$KONT_TOTAL <- rowSums(grades_no_na[, Mi_or_Zi_or_Lab])
boxplot(KONT_TOTAL~Grupa, grades_no_na, main="Ukupni brojevi bodova kontinuirane nastave", col="azure")
means <- tapply(grades_no_na$KONT_TOTAL,grades_no_na$Grupa,mean)
points(means,col="coral1",pch=16)

grades_no_na %>%
  group_by(Grupa)%>%
  summarise(mean = mean(KONT_TOTAL), median=median(KONT_TOTAL), std=sd(KONT_TOTAL))

Iz box plota možemo zaključiti da grupa 1 postiže bolje rezultate od na kontinuiranoj nastavi od grupa 2 i 3. Narančasta točka predstavlja srednju vrijednost a crna crta medijan. Po tim mjerama grupe 2 i 3 su otprilike iste. Ipak desni rep grupe 2 jest teži od grupe 3.

grades_no_na_IR <- grades[!(rowSums(is.na(grades[ ,just_IR])) > 0), ]
grades_no_na_IR$IR_TOTAL <- rowSums(grades_no_na_IR[, just_IR])
boxplot(IR_TOTAL~Grupa, grades_no_na_IR, main="Ukupni brojevi bodova na ispitnom roku", col="azure")
means <- tapply(grades_no_na_IR$IR_TOTAL,grades_no_na_IR$Grupa,mean)
points(means,col="coral1",pch=16)

grades_no_na_IR %>%
  group_by(Grupa)%>%
  summarise(mean = mean(IR_TOTAL), median=median(IR_TOTAL), std=sd(IR_TOTAL))

Poredak uspješnosti grupa se obrće na ispitnom roku. Grupe 2 i 3 definitivno postižu bolje rezultate od grupe 1, no teško je reći je li prednost grupe 3 nad grupom 2 značajna.

• Postoji li povezanost između uspjeha studenata na međuispitu i završnom ispitu (vrijedi li da su uspješniji studenti na MI ujedno uspješniji i na ZI)?

MI_and_ZI <- colnames(grades)[colnames(grades) %like% 'M|ZI']
grades_MI_ZI <- grades[!(rowSums(is.na(grades[ ,MI_and_ZI])) > 0), ]
grades_MI_ZI$MI_TOTAL <- rowSums(grades_MI_ZI[, just_MI])
grades_MI_ZI$ZI_TOTAL <- rowSums(grades_MI_ZI[, just_ZI])
plot(grades_MI_ZI$MI_TOTAL, grades_MI_ZI$ZI_TOTAL, main='Bodovi na MI i ZI', xlab='MI', ylab='ZI', col='darkblue')

cor(grades_MI_ZI$MI_TOTAL, grades_MI_ZI$ZI_TOTAL)

## [1] 0.06734055

Iz scatter plota te niske korelacije ove dvije varijable možemo zaključiti da ne postoji povezanost između uspjeha ispitima MI i ZI.

• Postoji li povezanost između uspjeha studenata na nekim zadatcima na ispitima i pojedinim laboratorijskim vježbama? Razmislite koji su mogući uzroci ovakvih zavisnosti, ako postoje.

Da, već smo utvrdili da u zadatku korelacijska analiza postoje određene korelacije između uspjeha na nekim laboratorijskim vježbama i određenim ispitnim zadatcima. Točnije, uspjeh na LAB_1 koreliran je s uspjehom na MI_1, a uspjeh na LAB_2 s uspjehom na ZI_1 te je to uočeno na cijeloj studentskoj populaciji. Dodatno, među studentima koji su izašli na ispitni rok postoji negativna korelacija između LAB_1 i ZI_1. Obrazloženje pozitivnih korelacija jest preklapanje gradiva labosa i zadataka. Teže je utvrditi zašto u određenim slučajevima dolazi do negativne korelacije.

Vizualizirajmo neke od najjačih korelacija:

x <- na.omit(grades[, c('LAB_1', 'MI_1')])
x_cor = round(cor(x$LAB_1, x$MI_1), 3)
ggplot(x, aes(x=LAB_1, y=MI_1) ) +
  stat_density_2d(aes(fill = ..density..), geom = "raster", contour = FALSE) +
  scale_fill_distiller(palette= "Spectral", direction=-1) +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  ggtitle(glue::glue('Bodovi na LAB_1 i MI_1, cor={x_cor}')) +
  theme(plot.title = element_text(hjust = 0.5))

x <- na.omit(grades[, c('LAB_2', 'ZI_1')])
x_cor = round(cor(x$LAB_2, x$ZI_1), 3)
ggplot(x, aes(x=LAB_2, y=ZI_1) ) +
  stat_density_2d(aes(fill = ..density..), geom = "raster", contour = FALSE) +
  scale_fill_distiller(palette= "Spectral", direction=-1) +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  ggtitle(glue::glue('Bodovi na LAB_2 i ZI_1, cor={x_cor}')) +
  theme(plot.title = element_text(hjust = 0.5))

x <- na.omit(grades[, c('LAB_1', 'ZI_1', just_IR)])
x_cor = round(cor(x$LAB_1, x$ZI_1), 3)
ggplot(x, aes(x=LAB_1, y=ZI_1) ) +
  stat_density_2d(aes(fill = ..density..), geom = "raster", contour = FALSE) +
  scale_fill_distiller(palette= "Spectral", direction=-1) +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  ggtitle(glue::glue('Bodovi na LAB_1 i ZI_1 među studentima koji su izašli na IR, cor={x_cor}')) +
  theme(plot.title = element_text(hjust = 0.5))

Postavite i analizirajte na ovaj način još barem jedno vlastito istraživačko pitanje.

Prati li ukupan broj bodova na ZI, MI i laboratorijskim vježbama normalnu distribuciju?

x <- na.omit(grades[, Mi_or_Zi_or_Lab])
x$TOTAL <- rowSums(x[, Mi_or_Zi_or_Lab])
hist(x$TOTAL, col="azure")

qqnorm(x$TOTAL, pch = 1, frame = FALSE)
qqline(x$TOTAL, col = "steelblue", lwd = 2)

skewness(x$TOTAL)

## [1] -0.304563

kurtosis(x$TOTAL)

## [1] 3.024189

shapiro.test(x$TOTAL)

## 
##  Shapiro-Wilk normality test
## 
## data:  x$TOTAL
## W = 0.99287, p-value = 0.01915

ad.test(x$TOTAL)

## 
##  Anderson-Darling normality test
## 
## data:  x$TOTAL
## A = 0.85976, p-value = 0.0271

Na temelju histograma i qqplota možemo zaključiti da je distribucija ukupnog zbroja zakrivljena u lijevo. Kurtosis distribucije poprilično dobro odgovara normalnoj distribuciji. Iz rezultata Shapiro–Wilk i Anderson–Darling testa možemo odbaciti hipotezu da uzorak dolazi iz normalne distribucije s 0.05 razinom značajnosti no ne i s 0.01. Ipak iz qqplota vidimo da odstupanje nije pretjerano izraženo pa u daljnjoj analizi pretpostavka normalnosti ne bi trebala biti previše narušena.

Prate li brojevi bodova LAB_2, i LAB_1 multivarijatnu normalnu distribuciju?

x <- na.omit(grades[, c('LAB_1', 'LAB_2')])
ggplot(x, aes(x=x$LAB_1, y=x$LAB_2))+
  geom_point(alpha = .2) +
  geom_density_2d()+
  theme_bw()

cov(x)

##          LAB_1    LAB_2
## LAB_1 1.001000 0.327244
## LAB_2 0.327244 1.009142

mvnorm.etest(x, R=200)

## 
##  Energy test of multivariate normality: estimated parameters
## 
## data:  x, sample size 493, dimension 2, replicates 200
## E-statistic = 1.7343, p-value < 2.2e-16

dens <- kde2d(x$LAB_1, x$LAB_2)

fig <- plot_ly(x = dens$x,
        y = dens$y,
        z = dens$z) %>% add_surface()
fig

Bodovi očito nisu normalno distribuirani s obzirom da poprimaju diskretne vrijednosti, te to potvrđuje i niska p vrijednost E-testa.

4.2. Regresijska analiza

Razmotrimo u kakvom su odnosu zadatci ispitnog roka s ostalim aktivnostima iz kontinuirane nastave. Istražite odnos koristeći model multivarijatne linearne regresije. Procijenite model gdje su zavisne varijable bodovi zadataka s ispitnog roka, odaberite konačni skup ulaznih varijabli i provjerite adekvatnost modela.

all_except_Grupa <- colnames(grades)[!(colnames(grades) %like% 'G')]
x <- na.omit(grades[, all_except_Grupa])
mlm.fit <- lm(cbind(IR_1, IR_2, IR_3, IR_4, IR_5) ~ ., data = x)
summary(mlm.fit)

## Response IR_1 :
## 
## Call:
## lm(formula = IR_1 ~ MI_1 + MI_2 + MI_3 + MI_4 + MI_5 + LAB_1 + 
##     ZI_1 + ZI_2 + ZI_3 + ZI_4 + ZI_5 + LAB_2, data = x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -17.8234  -1.2973  -0.1393   1.4238   5.2778 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  9.766210   5.025698   1.943   0.0555 .  
## MI_1         2.117073   0.483963   4.374 3.59e-05 ***
## MI_2        -0.490113   0.215913  -2.270   0.0259 *  
## MI_3         0.218906   0.188025   1.164   0.2477    
## MI_4         0.284559   0.365805   0.778   0.4389    
## MI_5        -0.264495   0.224464  -1.178   0.2421    
## LAB_1       -0.768253   0.371780  -2.066   0.0420 *  
## ZI_1        -0.149721   0.260230  -0.575   0.5667    
## ZI_2        -0.087315   0.395259  -0.221   0.8257    
## ZI_3        -0.008015   0.219966  -0.036   0.9710    
## ZI_4        -0.140767   0.353895  -0.398   0.6919    
## ZI_5         0.078622   0.354919   0.222   0.8252    
## LAB_2       -0.135750   0.433779  -0.313   0.7551    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.024 on 81 degrees of freedom
## Multiple R-squared:  0.2674, Adjusted R-squared:  0.1589 
## F-statistic: 2.464 on 12 and 81 DF,  p-value: 0.008531
## 
## 
## Response IR_2 :
## 
## Call:
## lm(formula = IR_2 ~ MI_1 + MI_2 + MI_3 + MI_4 + MI_5 + LAB_1 + 
##     ZI_1 + ZI_2 + ZI_3 + ZI_4 + ZI_5 + LAB_2, data = x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11.2134  -1.5310   0.2128   1.7150   5.0261 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  8.01148    4.59808   1.742 0.085242 .  
## MI_1        -1.05180    0.44278  -2.375 0.019891 *  
## MI_2         0.21372    0.19754   1.082 0.282514    
## MI_3        -0.04998    0.17203  -0.291 0.772138    
## MI_4         1.27026    0.33468   3.795 0.000283 ***
## MI_5         0.40002    0.20536   1.948 0.054896 .  
## LAB_1        0.74989    0.34015   2.205 0.030317 *  
## ZI_1         0.24192    0.23809   1.016 0.312600    
## ZI_2         0.20163    0.36163   0.558 0.578687    
## ZI_3        -0.15243    0.20125  -0.757 0.450985    
## ZI_4        -0.28305    0.32378  -0.874 0.384588    
## ZI_5         0.19543    0.32472   0.602 0.548961    
## LAB_2        0.17482    0.39687   0.440 0.660759    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.767 on 81 degrees of freedom
## Multiple R-squared:  0.258,  Adjusted R-squared:  0.1481 
## F-statistic: 2.347 on 12 and 81 DF,  p-value: 0.01221
## 
## 
## Response IR_3 :
## 
## Call:
## lm(formula = IR_3 ~ MI_1 + MI_2 + MI_3 + MI_4 + MI_5 + LAB_1 + 
##     ZI_1 + ZI_2 + ZI_3 + ZI_4 + ZI_5 + LAB_2, data = x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11.3259  -0.8119   0.1808   0.8788   3.4027 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  9.47499    3.15698   3.001  0.00357 ** 
## MI_1        -0.03460    0.30401  -0.114  0.90967    
## MI_2         0.11666    0.13563   0.860  0.39225    
## MI_3        -0.37164    0.11811  -3.147  0.00231 ** 
## MI_4         0.45839    0.22979   1.995  0.04942 *  
## MI_5         0.11588    0.14100   0.822  0.41358    
## LAB_1        1.00737    0.23354   4.313 4.50e-05 ***
## ZI_1         0.92339    0.16347   5.649 2.33e-07 ***
## ZI_2        -0.35880    0.24829  -1.445  0.15229    
## ZI_3        -0.16664    0.13818  -1.206  0.23133    
## ZI_4        -0.45180    0.22231  -2.032  0.04540 *  
## ZI_5        -0.06771    0.22295  -0.304  0.76214    
## LAB_2       -1.16383    0.27249  -4.271 5.26e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.9 on 81 degrees of freedom
## Multiple R-squared:  0.4181, Adjusted R-squared:  0.3319 
## F-statistic:  4.85 on 12 and 81 DF,  p-value: 5.902e-06
## 
## 
## Response IR_4 :
## 
## Call:
## lm(formula = IR_4 ~ MI_1 + MI_2 + MI_3 + MI_4 + MI_5 + LAB_1 + 
##     ZI_1 + ZI_2 + ZI_3 + ZI_4 + ZI_5 + LAB_2, data = x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9.8038 -2.8049  0.0276  2.5489 10.7682 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -5.24754    7.04645  -0.745   0.4586    
## MI_1        -0.03872    0.67856  -0.057   0.9546    
## MI_2        -0.31807    0.30273  -1.051   0.2965    
## MI_3         0.19983    0.26363   0.758   0.4506    
## MI_4         0.88201    0.51289   1.720   0.0893 .  
## MI_5        -0.08089    0.31472  -0.257   0.7978    
## LAB_1        0.48543    0.52127   0.931   0.3545    
## ZI_1        -0.34393    0.36486  -0.943   0.3487    
## ZI_2         2.43480    0.55419   4.393 3.35e-05 ***
## ZI_3         0.12054    0.30841   0.391   0.6969    
## ZI_4        -1.01527    0.49619  -2.046   0.0440 *  
## ZI_5        -0.27154    0.49763  -0.546   0.5868    
## LAB_2        0.94886    0.60819   1.560   0.1226    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.24 on 81 degrees of freedom
## Multiple R-squared:  0.2925, Adjusted R-squared:  0.1877 
## F-statistic:  2.79 on 12 and 81 DF,  p-value: 0.003117
## 
## 
## Response IR_5 :
## 
## Call:
## lm(formula = IR_5 ~ MI_1 + MI_2 + MI_3 + MI_4 + MI_5 + LAB_1 + 
##     ZI_1 + ZI_2 + ZI_3 + ZI_4 + ZI_5 + LAB_2, data = x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.8061 -0.7228  0.0193  1.0524  3.9388 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.127228   3.005118   0.708   0.4811    
## MI_1        -0.134391   0.289386  -0.464   0.6436    
## MI_2         0.191124   0.129105   1.480   0.1427    
## MI_3         0.015791   0.112430   0.140   0.8886    
## MI_4         0.164076   0.218733   0.750   0.4554    
## MI_5         0.112145   0.134218   0.836   0.4059    
## LAB_1        0.005849   0.222306   0.026   0.9791    
## ZI_1        -0.030631   0.155605  -0.197   0.8444    
## ZI_2         0.191384   0.236346   0.810   0.4204    
## ZI_3         0.029829   0.131529   0.227   0.8212    
## ZI_4         0.361210   0.211612   1.707   0.0917 .  
## ZI_5         1.208432   0.212224   5.694 1.92e-07 ***
## LAB_2       -0.290503   0.259378  -1.120   0.2660    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.808 on 81 degrees of freedom
## Multiple R-squared:  0.3326, Adjusted R-squared:  0.2338 
## F-statistic: 3.364 on 12 and 81 DF,  p-value: 0.000526

S obzirom na izrazito niske vrijednosti R-squared mjere možemo zaključiti da naši modeli loše objašnjavaju varijabilnost u podatcima, najviše oko 40%. Najjače veze u podatcima pronalazimo tako da pogledamo regresore s najvećim vrijednostima koeficijenata:

IR_1 ~ MI_1
IR_2 ~ MI_4
IR_3 ~ LAB_1, ZI_1, ↓LAB_2
IR_4 ~ ZI_2
IR_5 ~ ZI_5

Ove veze u skladu su s rezultatima korelacijske analize u kojoj smo utvrdili da su dotične varijable značajno korelirane.

Provjerimo imaju li stvarno svi regresori značajan utjecaj na model na razini signifikantnosti od 1%:

Anova(mlm.fit)

## 
## Type II MANOVA Tests: Pillai test statistic
##       Df test stat approx F num Df den Df    Pr(>F)    
## MI_1   1   0.53759  17.9039      5     77 9.741e-12 ***
## MI_2   1   0.33363   7.7105      5     77 6.392e-06 ***
## MI_3   1   0.41646  10.9908      5     77 5.263e-08 ***
## MI_4   1   0.24011   4.8662      5     77 0.0006395 ***
## MI_5   1   0.23746   4.7955      5     77 0.0007206 ***
## LAB_1  1   0.59924  23.0267      5     77 4.612e-14 ***
## ZI_1   1   0.62610  25.7873      5     77 3.398e-15 ***
## ZI_2   1   0.39508  10.0580      5     77 1.953e-07 ***
## ZI_3   1   0.09760   1.6657      5     77 0.1529827    
## ZI_4   1   0.14903   2.6971      5     77 0.0267456 *  
## ZI_5   1   0.32774   7.5078      5     77 8.754e-06 ***
## LAB_2  1   0.42316  11.2973      5     77 3.451e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

P-vrijednost testa je dovoljno velika da ne možemo odbaciti nultu hipotezu da kompleksniji model ne donosi značajno poboljšanje u odnosu na jednostavniji model (raz. sign. od 1%). Prema tome, odlučujemo se za model bez regresora ZI_3, te ZI_4.

Izbacimo ZI_3, i ZI_4 te ponovimo analizu:

mlm2.fit <- lm(cbind(IR_1, IR_2, IR_3, IR_4, IR_5) ~ . - ZI_3 - ZI_4, data = x)
anova(mlm.fit, mlm2.fit)

lh.out <- linearHypothesis(mlm.fit, hypothesis.matrix = c("ZI_3 = 0", "ZI_4 = 0"))
lh.out

## 
## Sum of squares and products for the hypothesis:
##           IR_1      IR_2       IR_3      IR_4       IR_5
## IR_1  1.447222  2.941923   4.677857  10.39632  -3.715770
## IR_2  2.941923  9.513261  13.161183  16.51156  -7.793521
## IR_3  4.677857 13.161183  18.895365  28.82607 -12.258660
## IR_4 10.396322 16.511558  28.826066  80.73071 -26.378665
## IR_5 -3.715770 -7.793521 -12.258660 -26.37867   9.556628
## 
## Sum of squares and products for error:
##          IR_1     IR_2      IR_3       IR_4      IR_5
## IR_1 740.7699 409.8653 307.68061  247.55419 117.06089
## IR_2 409.8653 620.0737 137.13962  278.37382 104.56566
## IR_3 307.6806 137.1396 292.30326  338.64132  93.83889
## IR_4 247.5542 278.3738 338.64132 1456.23410   6.63055
## IR_5 117.0609 104.5657  93.83889    6.63055 264.85834
## 
## Multivariate Tests: 
##                  Df test stat approx F num Df den Df    Pr(>F)   
## Pillai            2 0.2285414 2.012605     10    156 0.0353865 * 
## Wilks             2 0.7800254 2.036782     10    154 0.0330844 * 
## Hotelling-Lawley  2 0.2710267 2.059803     10    152 0.0310359 * 
## Roy               2 0.2214267 3.454257      5     78 0.0071522 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Temeljem p-vrijednosti provedenih testova dolazimo do istog zaključka, tj. da se trebamo odlučiti za jednostavniji model.

Provjerimo još značajnost koeficijenata jednostavnijeg modela:

Anova(mlm2.fit)

## 
## Type II MANOVA Tests: Pillai test statistic
##       Df test stat approx F num Df den Df    Pr(>F)    
## MI_1   1   0.50879  16.3652      5     79 4.699e-11 ***
## MI_2   1   0.29873   6.7307      5     79 2.844e-05 ***
## MI_3   1   0.34704   8.3975      5     79 2.089e-06 ***
## MI_4   1   0.21853   4.4184      5     79  0.001339 ** 
## MI_5   1   0.21800   4.4045      5     79  0.001372 ** 
## LAB_1  1   0.55167  19.4418      5     79 1.430e-12 ***
## ZI_1   1   0.58149  21.9526      5     79 1.017e-13 ***
## ZI_2   1   0.44491  12.6636      5     79 4.857e-09 ***
## ZI_5   1   0.28926   6.4302      5     79 4.624e-05 ***
## LAB_2  1   0.39650  10.3805      5     79 1.122e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Svi regresori u jednostavnijem modelu značajni su na razini signifikantnosti 1%.

Kako bismo provjerili adekvatnost modela proučit ćemo grafove reziduala:

## Preuzeto iz auditorne vjezbe 4
## define our own "rstandard" method for "mlm" class
rstandard.mlm2 <- function (model) {
  Q <- with(model, qr.qy(qr, diag(1, nrow = nrow(qr$qr), ncol = qr$rank)))  ## Q matrix
  hii <- rowSums(Q ^ 2)  ## diagonal of hat matrix QQ'
  RSS <- colSums(model$residuals ^ 2)  ## residual sums of squares (for each model)
  sigma <- sqrt(RSS / model$df.residual)  ##  ## Pearson estimate of residuals (for each model)
  pointwise_sd <- outer(sqrt(1 - hii), sigma)  ## point-wise residual standard error (for each model)
  model$residuals / pointwise_sd  ## standardised residuals
  }


plot(fitted(mlm2.fit), rstandard(mlm2.fit), col = as.numeric(col(fitted(mlm2.fit))), pch = 19)
abline(h=0)
legend("bottomright", legend = paste0("IR ", 1:ncol(fitted(mlm2.fit))), pch = 19,
       col = 1:ncol(fitted(mlm2.fit)), text.col = 1:ncol(fitted(mlm2.fit)))

Generalno reziduali izgledaju nasumično raspršeni oko nule što ukazuje na dobar fit modela. Ipak, primjećujemo outlier u grupama reziduala IR_1, IR_2, te IR_3. Još jedno odstupanje od očekivanog je ljevkast oblik reziduala IR_1. Drugim riječima s porastom predviđene vrijednosti IR_1 varijanca reziduala pada.

To možemo bolje vidjeti na sljedećem grafu:

plot(fitted(mlm2.fit)[, 'IR_1'], rstandard(mlm2.fit)[, 'IR_1'], col = 'darkblue', pch = 19)
abline(h=0)

Pogledajmo kojim retcima iz matrice grades pripadaju outlieri reziduala:

st_r <- rstandard(mlm2.fit)
which(st_r[, c('IR_2', 'IR_3', 'IR_1')] < -3, arr.ind = TRUE)

##    row col
## 35   8   1
## 35   8   2
## 35   8   3

grades[35,]

Vidimo da su svi outlieri potekli od istog retka kojeg smo već prije detektirali kao outlier. U pitanju je redak koji ima vrijednost 0 za sve zadatke na ispitnom roku. Iako ovaj redak vjerojatno nije pogrešan izbacit ćemo ga iz našeg modela jer ne želimo da previše “povuče” koeficijente regresora.

grades2 <- grades[-35,]
x2 <- na.omit(grades2[, all_except_Grupa])
mlm_f.fit <- lm(cbind(IR_1, IR_2, IR_3, IR_4, IR_5) ~ . - ZI_3 - ZI_4, data = x2)
summary(mlm_f.fit)

## Response IR_1 :
## 
## Call:
## lm(formula = IR_1 ~ (MI_1 + MI_2 + MI_3 + MI_4 + MI_5 + LAB_1 + 
##     ZI_1 + ZI_2 + ZI_3 + ZI_4 + ZI_5 + LAB_2) - ZI_3 - ZI_4, 
##     data = x2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.2544 -1.3273 -0.1423  1.1450  4.8465 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 13.26137    3.60192   3.682 0.000414 ***
## MI_1         2.50607    0.33885   7.396 1.08e-10 ***
## MI_2        -0.66385    0.15498  -4.283 4.97e-05 ***
## MI_3         0.22849    0.12941   1.766 0.081184 .  
## MI_4         0.03733    0.26464   0.141 0.888160    
## MI_5        -0.39847    0.16216  -2.457 0.016106 *  
## LAB_1       -1.19467    0.27154  -4.400 3.24e-05 ***
## ZI_1        -0.40045    0.18423  -2.174 0.032621 *  
## ZI_2        -0.28981    0.27063  -1.071 0.287369    
## ZI_5         0.08105    0.25276   0.321 0.749294    
## LAB_2        0.18514    0.31468   0.588 0.557914    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.183 on 82 degrees of freedom
## Multiple R-squared:  0.4921, Adjusted R-squared:  0.4302 
## F-statistic: 7.945 on 10 and 82 DF,  p-value: 8.355e-09
## 
## 
## Response IR_2 :
## 
## Call:
## lm(formula = IR_2 ~ (MI_1 + MI_2 + MI_3 + MI_4 + MI_5 + LAB_1 + 
##     ZI_1 + ZI_2 + ZI_3 + ZI_4 + ZI_5 + LAB_2) - ZI_3 - ZI_4, 
##     data = x2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.7744 -1.5268 -0.0431  1.6066  4.4480 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  9.596262   3.991002   2.404 0.018449 *  
## MI_1        -0.707022   0.375455  -1.883 0.063231 .  
## MI_2         0.079412   0.171721   0.462 0.644984    
## MI_3        -0.002924   0.143392  -0.020 0.983780    
## MI_4         1.095121   0.293226   3.735 0.000346 ***
## MI_5         0.309741   0.179674   1.724 0.088493 .  
## LAB_1        0.444296   0.300873   1.477 0.143589    
## ZI_1         0.041729   0.204135   0.204 0.838531    
## ZI_2        -0.028507   0.299861  -0.095 0.924492    
## ZI_5         0.244872   0.280062   0.874 0.384481    
## LAB_2        0.394770   0.348674   1.132 0.260848    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.418 on 82 degrees of freedom
## Multiple R-squared:  0.2408, Adjusted R-squared:  0.1483 
## F-statistic: 2.601 on 10 and 82 DF,  p-value: 0.008559
## 
## 
## Response IR_3 :
## 
## Call:
## lm(formula = IR_3 ~ (MI_1 + MI_2 + MI_3 + MI_4 + MI_5 + LAB_1 + 
##     ZI_1 + ZI_2 + ZI_3 + ZI_4 + ZI_5 + LAB_2) - ZI_3 - ZI_4, 
##     data = x2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.7007 -0.7963  0.0838  0.7797  3.1684 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 10.859458   2.255468   4.815 6.64e-06 ***
## MI_1         0.327232   0.212184   1.542 0.126874    
## MI_2        -0.020053   0.097046  -0.207 0.836805    
## MI_3        -0.317817   0.081036  -3.922 0.000182 ***
## MI_4         0.282825   0.165713   1.707 0.091660 .  
## MI_5         0.028964   0.101540   0.285 0.776178    
## LAB_1        0.680437   0.170035   4.002 0.000137 ***
## ZI_1         0.720743   0.115364   6.248 1.76e-08 ***
## ZI_2        -0.647450   0.169463  -3.821 0.000258 ***
## ZI_5         0.001458   0.158274   0.009 0.992673    
## LAB_2       -0.927088   0.197049  -4.705 1.02e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.367 on 82 degrees of freedom
## Multiple R-squared:  0.4737, Adjusted R-squared:  0.4095 
## F-statistic: 7.381 on 10 and 82 DF,  p-value: 3.108e-08
## 
## 
## Response IR_4 :
## 
## Call:
## lm(formula = IR_4 ~ (MI_1 + MI_2 + MI_3 + MI_4 + MI_5 + LAB_1 + 
##     ZI_1 + ZI_2 + ZI_3 + ZI_4 + ZI_5 + LAB_2) - ZI_3 - ZI_4, 
##     data = x2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.3106 -3.1140 -0.0312  2.9942 10.7700 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -4.0121     6.8500  -0.586 0.559673    
## MI_1          0.1273     0.6444   0.198 0.843876    
## MI_2         -0.3844     0.2947  -1.304 0.195787    
## MI_3          0.1908     0.2461   0.775 0.440473    
## MI_4          0.7850     0.5033   1.560 0.122687    
## MI_5         -0.1141     0.3084  -0.370 0.712262    
## LAB_1         0.1712     0.5164   0.332 0.741075    
## ZI_1         -0.4221     0.3504  -1.205 0.231720    
## ZI_2          2.0604     0.5147   4.003 0.000136 ***
## ZI_5         -0.1992     0.4807  -0.414 0.679718    
## LAB_2         1.2063     0.5984   2.016 0.047098 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.151 on 82 degrees of freedom
## Multiple R-squared:   0.27,  Adjusted R-squared:  0.181 
## F-statistic: 3.033 on 10 and 82 DF,  p-value: 0.002585
## 
## 
## Response IR_5 :
## 
## Call:
## lm(formula = IR_5 ~ (MI_1 + MI_2 + MI_3 + MI_4 + MI_5 + LAB_1 + 
##     ZI_1 + ZI_2 + ZI_3 + ZI_4 + ZI_5 + LAB_2) - ZI_3 - ZI_4, 
##     data = x2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.0875 -0.7469 -0.0542  1.0434  4.0554 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.799946   2.840375   1.338    0.185    
## MI_1        -0.041869   0.267210  -0.157    0.876    
## MI_2         0.137374   0.122213   1.124    0.264    
## MI_3         0.003128   0.102051   0.031    0.976    
## MI_4         0.080482   0.208687   0.386    0.701    
## MI_5         0.058133   0.127873   0.455    0.651    
## LAB_1       -0.088081   0.214130  -0.411    0.682    
## ZI_1        -0.110793   0.145282  -0.763    0.448    
## ZI_2         0.263702   0.213409   1.236    0.220    
## ZI_5         1.161787   0.199319   5.829 1.06e-07 ***
## LAB_2       -0.221749   0.248149  -0.894    0.374    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.721 on 82 degrees of freedom
## Multiple R-squared:   0.32,  Adjusted R-squared:  0.237 
## F-statistic: 3.858 on 10 and 82 DF,  p-value: 0.0002624

plot(fitted(mlm_f.fit), rstandard(mlm_f.fit), col = as.numeric(col(fitted(mlm_f.fit))), pch = 19)
abline(h=0)
legend("bottomright", legend = paste0("IR ", 1:ncol(fitted(mlm_f.fit))), pch = 19,
       col = 1:ncol(fitted(mlm_f.fit)), text.col = 1:ncol(fitted(mlm_f.fit)))

Standardizirani reziduali još uvijek izgledaju kao nasumičan pojas oko nule, osim već zamijećene deformacije IR_1 reziduala.

Provjerimo normalnost reziduala:

st_r <- rstandard(mlm_f.fit)
for (val in colnames(st_r)) {
  qqnorm(st_r[,val], pch = 1, frame = FALSE, main=glue::glue('{val} residuals QQ plot'))
  qqline(st_r[,val], col = "steelblue", lwd = 2)
}

Iz QQ plotova možemo zaključiti da su pretpostavke o normalnoj distribuiranosti reziduala u određenoj mjeri prekršene:

IR_1 -> pretežak desni rep
IR_2 -> dobro slaganje s teoretskom distribucijom
IR_3 -> malo pretežak lijevi rep
IR_4 -> malo pretanki repovi
IR_5 -> pretežak lijevi rep

Iz toga možemo zaključiti da ne možemo računati intervali pouzdanosti bez određene mjere opreza. Moguće je da, zbog loših temelja, ne reflektiraju stvarnu situaciju.

Ovakav rezultat nije iznenađujuć s obzirom da naš model objašnjava relativno malo varijance podataka. Trebali bi potražiti još relevantnih regresora. Npr prosjek, uspjeh na ostalim kolegijima, grupa itd.

Najgore odstupanje od normalnosti primjećujemo za IR_1 i IR_5 reziduale stoga moramo biti najoprezniji u predikciji tih varijabli.

Provjerimo zavisnost reziduala IR_5 i IR_1 o nezavisnim varijablama i prediktivnoj vrijednosti. IR_5:

just_MI <- colnames(grades)[colnames(grades) %like% 'MI']
just_rel_ZI <- c('ZI_1', 'ZI_2', 'ZI_5')
just_IR <- colnames(grades)[colnames(grades) %like% 'IR']
just_LAB <- colnames(grades)[colnames(grades) %like% 'LAB']
par(mfrow=c(3,2), mar=c(3.3,3.1,1,0), mgp=c(1.5, 0.5, 0))
for (val in just_MI) {
  plot(x2[, val], rstandard(mlm_f.fit)[, 'IR_5'], xlab=val, ylab='Standardized residuals', col = 'darkblue', pch = 19, main=glue::glue('IR_5 residuals by {val}'))
  abline(h=0)
}

par(mfrow=c(3,2), mar=c(3.3,3.1,1,0), mgp=c(1.5, 0.5, 0))

for (val in just_rel_ZI) {
  plot(x2[, val], rstandard(mlm_f.fit)[, 'IR_5'], xlab=val, ylab='Standardized residuals', col = 'darkblue', pch = 19, main=glue::glue('IR_5 residuals by {val}'))
  abline(h=0)
}

par(mfrow=c(2,1), mar=c(3.3,3.1,1,0), mgp=c(1.5, 0.5, 0))

for (val in just_LAB) {
  plot(x2[, val], rstandard(mlm_f.fit)[, 'IR_5'], xlab=val, ylab='Standardized residuals', col = 'darkblue', pch = 19, main=glue::glue('IR_5 residuals by {val}'))
  abline(h=0)
}

par(mfrow=c(3,2), mar=c(3.3,3.1,1,0), mgp=c(1.5, 0.5, 0))
for (val in colnames(fitted(mlm_f.fit))) {
  plot(fitted(mlm_f.fit)[, val], rstandard(mlm_f.fit)[, 'IR_5'], xlab=val, ylab='Standardized residuals', col = 'darkblue', pch = 19, main=glue::glue('IR_5 residuals by predicted {val}'))
  abline(h=0)
}

IR_1:

par(mfrow=c(3,2), mar=c(3.3,3.1,1,0), mgp=c(1.5, 0.5, 0))
for (val in just_MI) {
  plot(x2[, val], rstandard(mlm_f.fit)[, 'IR_1'], xlab=val, ylab='Standardized residuals', col = 'darkblue', pch = 19, main=glue::glue('IR_1 residuals by {val}'))
  abline(h=0)
}

par(mfrow=c(3,2), mar=c(3.3,3.1,1,0), mgp=c(1.5, 0.5, 0))

for (val in just_rel_ZI) {
  plot(x2[, val], rstandard(mlm_f.fit)[, 'IR_1'], xlab=val, ylab='Standardized residuals', col = 'darkblue', pch = 19, main=glue::glue('IR_1 residuals by {val}'))
  abline(h=0)
}

par(mfrow=c(2,1), mar=c(3.3,3.1,1,0), mgp=c(1.5, 0.5, 0))

for (val in just_LAB) {
  plot(x2[, val], rstandard(mlm_f.fit)[, 'IR_1'], xlab=val, ylab='Standardized residuals', col = 'darkblue', pch = 19, main=glue::glue('IR_1 residuals by {val}'))
  abline(h=0)
}

par(mfrow=c(3,2), mar=c(3.3,3.1,1,0), mgp=c(1.5, 0.5, 0))
for (val in colnames(fitted(mlm_f.fit))) {
  plot(fitted(mlm_f.fit)[, val], rstandard(mlm_f.fit)[, 'IR_1'], xlab=val, ylab='Standardized residuals', col = 'darkblue', pch = 19, main=glue::glue('IR_1 residuals by predicted {val}'))
  abline(h=0)
}

Ne primjećujemo nikakve problematične zavisnosti osim već uočene nekonstantnosti varijance reziduala IR_1 s obzirom na ^IR_1.

Provjerimo potječu li problematični reziduali IR_1 i IR_5 iz istog redaka (jer ih je otprilike jednako):

which(rstandard(mlm_f.fit)[, 'IR_5'] < -1, arr.ind = FALSE)

##  26  36  50  63  97 124 131 155 187 193 307 337 414 431 
##   6   8   9  14  18  23  25  29  36  39  59  69  79  83

which(rstandard(mlm_f.fit)[, 'IR_1'] > 1.1, arr.ind=FALSE)

##   4   5  50  62 112 113 131 213 359 395 413 447 481 483 500 
##   2   3   9  13  21  22  25  42  70  76  78  86  90  91  93

Ne pripadaju.

Naš model možemo dalje koristiti za predikcije:

newdata <- data.frame(MI_1=6,  MI_2=7,  MI_3=5,  MI_4=8,  MI_5=8,  LAB_1=9, ZI_1=2,  ZI_2=2,  ZI_3=3,  ZI_4=4,  ZI_5=1,  LAB_2=5)
predict(mlm_f.fit, newdata, interval="confidence") 

##       IR_1     IR_2     IR_3     IR_4     IR_5
## 1 10.77838 23.37811 15.22424 11.03109 5.201043

predict(mlm_f.fit, newdata, interval="predict")

##       IR_1     IR_2     IR_3     IR_4     IR_5
## 1 10.77838 23.37811 15.22424 11.03109 5.201043

lm_ir5.fit <- lm(IR_5 ~ MI_1 + MI_2 + MI_3 + MI_4 + MI_5 + LAB_1 + ZI_1 + ZI_2 + ZI_3 + ZI_4 + ZI_5 + LAB_2 - ZI_3 - ZI_4, data = x2)
summary(lm_ir5.fit)

## 
## Call:
## lm(formula = IR_5 ~ MI_1 + MI_2 + MI_3 + MI_4 + MI_5 + LAB_1 + 
##     ZI_1 + ZI_2 + ZI_3 + ZI_4 + ZI_5 + LAB_2 - ZI_3 - ZI_4, data = x2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.0875 -0.7469 -0.0542  1.0434  4.0554 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.799946   2.840375   1.338    0.185    
## MI_1        -0.041869   0.267210  -0.157    0.876    
## MI_2         0.137374   0.122213   1.124    0.264    
## MI_3         0.003128   0.102051   0.031    0.976    
## MI_4         0.080482   0.208687   0.386    0.701    
## MI_5         0.058133   0.127873   0.455    0.651    
## LAB_1       -0.088081   0.214130  -0.411    0.682    
## ZI_1        -0.110793   0.145282  -0.763    0.448    
## ZI_2         0.263702   0.213409   1.236    0.220    
## ZI_5         1.161787   0.199319   5.829 1.06e-07 ***
## LAB_2       -0.221749   0.248149  -0.894    0.374    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.721 on 82 degrees of freedom
## Multiple R-squared:   0.32,  Adjusted R-squared:  0.237 
## F-statistic: 3.858 on 10 and 82 DF,  p-value: 0.0002624

predict(lm_ir5.fit, newdata, interval="confidence") 

##        fit     lwr      upr
## 1 5.201043 2.07993 8.322157

predict(lm_ir5.fit, newdata, interval="predict")

##        fit      lwr      upr
## 1 5.201043 0.567941 9.834146

Nažalost u R-u još ne postoje gotove funkcije za računanje confidence elipsoide oko predikcija, no s obzirom da naš model objašnjava malo od ukupne varijance podataka možemo očekivati jako veliku elipsoidu pouzdanosti. To u praksi znači da nam model nije koristan u predviđanju budućih vrijednosti. Za ilustraciju izračunali smo interval pouzdanosti sredine i interval pouzdanosti predikcije samo za zavisnu varijablu IR_5 za neki naš izmišljen primjer. Vidimo da nam ti intervali nisu uopće korisni jer su u oba slučaja značajno širi od interkvartilnog raspona distribucije IR_5 [5.0, 7.5].